The response of a metapopulation to a changing environment

A species distributed across diverse environments may adapt to local conditions. We ask how quickly such a species changes its range in response to changed conditions. Szép et al. (Szép E, Sachdeva H, Barton NH. 2021 Polygenic local adaptation in metapopulations: a stochastic eco-evolutionary model. Evolution 75, 1030–1045 (doi:10.1111/evo.14210)) used the infinite island model to find the stationary distribution of allele frequencies and deme sizes. We extend this to find how a metapopulation responds to changes in carrying capacity, selection strength, or migration rate when deme sizes are fixed. We further develop a ‘fixed-state’ approximation. Under this approximation, polymorphism is only possible for a narrow range of habitat proportions when selection is weak compared to drift, but for a much wider range otherwise. When rates of selection or migration relative to drift change in a single deme of the metapopulation, the population takes a time of order m−1 to reach the new equilibrium. However, even with many loci, there can be substantial fluctuations in net adaptation, because at each locus, alleles randomly get lost or fixed. Thus, in a finite metapopulation, variation may gradually be lost by chance, even if it would persist in an infinite metapopulation. When conditions change across the whole metapopulation, there can be rapid change, which is predicted well by the fixed-state approximation. This work helps towards an understanding of how metapopulations extend their range across diverse environments. This article is part of the theme issue ‘Species’ ranges in the face of changing environments (Part II)’.


Introduction
Species must adapt to varied environments, while drawing on a common pool of genetic variation. Thus, there is a tension between selection that favours different alleles in different places, and the maintenance of diversity across the whole species. Local populations can only sustain themselves if they are sufficiently well-adapted; conversely, adaptation to conditions beyond the current niche can extend the range of the species.
These issues, which lie at the interface between ecology and evolution, have only quite recently attracted sustained theoretical attention. This ranges from studies of 'evolutionary rescue', typically of a single isolated deme [1][2][3], through to analyses of limits to a species' range in one or two spatial dimensions [4][5][6]. Here, we consider an idealized metapopulation; in this island model, there is no explicit spatial structure. Nevertheless, we can ask whether the species' range can extend over a variety of habitats, and examine how it responds dynamically to changing conditions-either in a single deme, or across the whole metapopulation.
We examine a simple model, in which directional selection favours alternative alleles in two different habitats. Provided that selection is stronger than migration, these alternative adaptations can be maintained despite gene flow. There is substantial literature on how heterogeneous selection can maintain diversity, beginning with Levene [7]. However, this is largely deterministic, neglecting random drift in small local populations. Here, we are primarily concerned with the erosion of adaptation by random drift within local demes-which can cause a substantial 'drift load' even when the whole metapopulation is very large. This paper is an extension of Szép et al. [8], which analysed the joint evolution of allele frequencies and deme sizes in an island model with explicit density-dependent regulation; a diffusion approximation gave explicit formulae for the stationary distribution of an infinite metapopulation. Here, we extend this treatment to consider the evolution of individual demes, and of the whole metapopulation, as conditions change (directly through changes in selection and gene flow, and indirectly through their effect on the population size); we also consider fluctuations in a metapopulation with a limited number of demes, where variation can be lost by chance. We simplify the problem by assuming that deme sizes are fixed, independent of adaptedness (soft selection), but believe that the methods we introduce can be extended to allow density regulation (hard selection).
In principle, we can calculate the joint distribution of deme size and allele frequencies under the diffusion approximation. However, this is numerically challenging, since it involves a high dimensional partial differential equation; in any case, it can only be done for an infinite metapopulation, where the mean population size and allele frequencies across the population as a whole change deterministically, even though population sizes and allele frequencies within any deme follow a distribution. In order to go beyond mere simulation, we use the approximation that loci are typically near fixation; this is accurate if the number of non-native alleles that enter per generation is small. It allows us to follow the distribution of states of a finite metapopulation through time, which depends only on the rates of substitutions in either direction. This 'fixed-state' approximation is an extension of models of 'adaptive walks' (e.g. [9,10]) to structured populations.
We first consider an infinite metapopulation, and determine the accuracy of the fixed-state approximation. We then apply the approximation to calculate the dynamics of a finite metapopulation, and to find how its equilibria depend on the number of demes. (In order for a non-trivial equilibrium to exist, we must allow a low rate of mutation to maintain variation in the long term.) Finally, we show how metapopulations respond to changing conditions, focusing on changes that take the system between qualitatively different regimes.

Model and methods
We simulate a haploid population, assuming linkage equilibrium. Provided that selection is weak, this is accurate, and allows us to efficiently simulate large numbers of loci and demes; Szép et al. [8,Supplementary Information C] examined the effects of linkage disequilibrium in this model, using individual-based simulations. We obtain analytical results by taking the diffusion limit, which also assumes weak selection, and then approximate this by assuming that demes are near fixation, which applies when there are few migrants (Nm < 1). As is traditional in population genetics, we take the fundamental model to be the diffusion, since this captures the behaviour of a variety of particular life histories, and identifies the key dimensionless parameters.

(a) Simulations
Our baseline island model assumes that demes each have carrying capacity of N haploid individuals, and contribute equally to the migrant pool. A deme of size N is expected to lose a fraction m of individuals by emigration, and receives a Poisson distributed number of migrants, Nm*, with expectation Nm. There are L biallelic loci, with the two alternative alleles labelled X i,k ¼ 0 or 1; i labels the deme, and k the locus. Deme i is described by {j i,1 , j i,2 , …, j i,L }, where 0 ≤ j i,k ≤ N is the number of copies of the '1' allele at the k 0 th locus. That allele is favoured by selection s i , which we assume to be the same across loci; the marginal relative fitnesses are 1 : e s i , and fitnesses multiply across loci. Under soft selection, loci evolve independently, and so it would be straightforward to extend to allow variation in selection across loci.
We assume linkage equilibrium (LE), and apply the Wright-Fisher model to each locus independently. After selection, allele frequencies are p Ã i,k ¼ j i,k =ððN À j i,k Þ e Às i þ j i,k Þ, and after migration, where p k is the frequency averaged across all demes of the metapopulation. The new population in deme i consists of N individuals, the number of allele copies at locus k being binomially sampled with frequency p ÃÃ i,k . This procedure is accurate provided that s is not too large (<0.2, say), so that recombination shuffles genes faster than selection, drift or migration build up associations between them [8, Supplementary Information C].
A Mathematica notebook containing the simulation code and result is provided as electronic supplementary material.

(b) Diffusion approximation
The diffusion approximation to this model describes the evolution of the joint distribution of allele frequencies across different demes, conditional on the mean allele frequency across the metapopulation [11]. A single deme follows a stochastic path governed by this distribution, while an infinite metapopulation represents the whole distribution, which evolves deterministically at the level of the metapopulation. The diffusion depends only on scaled parameters Ns, Nm.
Wright [12,13] gave an explicit solution for the stationary distribution of allele frequencies where Z is a normalizing constant and q = 1 − p. Under this simple model of directional selection, allele frequencies evolve independently across demes and across loci, conditional on the mean allele frequencies, p k . Equation (2.1) applies to a single deme; the subscript i was dropped for clarity. All demes that share the same parameters will follow the same distribution, in a given habitat, and so we can integrate over the distribution, and sum over habitats, to find the mean p k . This allows us to solve fully for the stationary state.

(i) Fixed-state approximation
If the number of incoming alleles is small (Nm ≪ 1) then the distribution of allele frequencies will be sharply peaked around 0 and 1. To a good approximation, populations are near fixation for one or other allele, and their state is determined by the rates of substitution in either direction. Since royalsocietypublishing.org/journal/rstb Phil. Trans. R. Soc. B 377: 20210009 we will later be considering the stationary state of a finite metapopulation, we must include mutation, which we assume to be symmetric at rate μ. Then, the rate at which demes currently fixed for allele 0 substitute allele 1, λ 0→1 (or vice versa, λ 1→0 ) is the product of the number of '1' (or '0') alleles entering the population, and their individual fixation probability (see electronic supplementary material for details). Thus: Different loci evolve independently, conditional on the numbers of migrants coming into the deme ðNm pÞ, ðNm qÞ. Note that when Ns and Nm are small, these rates are proportional to m in the neutral limit.
For an infinite metapopulation, and two habitats with selection s 1 , s 2 , with deme sizes fixed at N (i.e. soft selection), we can just follow the proportion of demes fixed for the '1' allele in each habitat, P 1 , P 2 . Neglecting mutation, the rates given by equation (2.2) lead immediately to @ t P 1 ¼ 2s 1 Nm 1 À e À2Ns1 ( pQ 1 À qe À2Ns1 P 1 ), @ t P 2 ¼ 2s 2 Nm 1 À e À2Ns2 ( pQ 2 À qe À2Ns2 P 2 ) and p ¼ r P 1 þ ð1 À rÞP 2 : ð2:3Þ The first two equations involve the difference in net rates of substitution in each direction, where Q, P are the fraction of loci near fixation for 0, 1; p, q are the fraction of migrants with allele 1 versus 0, which can contribute to a substitution; and the fixation probabilities in each direction are in the ratio 1 : e À2Ns1 . Finally, the mean allele frequencies are a weighted average across habitats, which are in the proportions ρ : 1 − ρ. These equations apply separately to each locus, but for simplicity, in numerical examples we will assume symmetric initial conditions, so that the proportion of demes fixed for the 1 allele in each habitat, P 1 , P 2 , are the same for all loci, and correspond to the proportion of loci fixed for the '1' allele in each deme.
If the '1' allele is favoured in habitat 1, but disfavoured in habitat 2 (i.e. s 2 < 0 < s 1 ), and if neither habitat is too rare, then polymorphism is possible, with equilibrium frequency given by p ¼ r(e 2N(s1Às2) À 1) À (e À2Ns2 À 1) (e 2Ns1 À 1)(e À2Ns2 À 1) , (as derived from equation (2.2); see electronic supplementary material). If selection is weak relative to drift, polymorphism is possible only for a very narrow range of habitat proportions (see left of figure 6 in appendix Aa), whereas if it is strong, polymorphism is possible over a wide range (right of figure 6 in appendix Aa). Suppose now that there are a finite number of demes, with d i having habitat i. At any one locus, the state of the metapopulation is described by the number of demes fixed for the '1' allele, 0 ≤ k i ≤ d i . For example, with two habitats, there are (d 1 + 1)(d 2 + 1) possible values for the state {k 1 , k 2 }. The probability of transitions between these states depends on the mean allele frequency across the metapopulation. With soft selection, where all demes have the same size N, this mean is just p ¼ ðk 1 þ k 2 Þ=ðd 1 þ d 2 Þ. We can therefore calculate the transition matrix that governs the stochastic evolution of the metapopulation; the stationary state is given by the leading eigenvector of this matrix. With soft selection, each locus evolves independently, governed by this matrix, and so we can easily calculate the stochastic evolution of the metapopulation.
In appendix Ab, we examine the accuracy of the fixedstate approximation under soft selection. This approximation applies in the limit of low migration, and identifies the failure of adaptation due to random drift.

Results (a) Evolution of a single deme
Consider a metapopulation, where Nm is small enough that populations are near fixation. If Ns 1 = 1 in a rare habitat, represented in ρ = 0.2 of the demes, and Ns 2 = −2 in the common habitat, then polymorphism will be maintained with p ¼ 0:079 overall (equation (2.4)). We begin by considering how a single deme responds to changes in its local conditions, for fixed p. In the focal deme, allele frequencies will be in the ratio q : p e 2Ns1 when Nm ≪ 1, since that is the ratio of substitution rates in either direction; hence, the expected allele frequency in the rare habitat is 0.386 (l.h.s. of figure 1a). As Nm increases, the expected allele frequency decreases, approaching p ¼ 0:079 (r.h.s. of figure 1a). For given Nm, the expected allele frequency in the focal deme increases with Ns 1 from p to 1, as selection becomes more effective (figure 1b). Figure 2a shows how the distribution of allele frequencies changes as Nm changes. If all loci start close to the frequency in the gene pool ð p ¼ 0:079Þ then with a low migration rate (Nm = 0.05), even weak selection (Ns = 1) can raise the mean substantially, to 0.355. However, this increase is slow, taking ∼5000 generations, because it occurs through occasional substitutions, at a rate proportional to m = 5 × 10 −4 (equation (2.3)). The population does mostly flip between fixation of one or other allele, giving a U-shaped frequency distribution (e.g. grey trajectory in figure 2), and so the fixed-state approximation is quite close to the exact mean (orange versus red at left). However, the average across even 100 loci fluctuates substantially (blue), implying that population fitness will fluctuate randomly, even when adaptation is highly polygenic.
At 20 000 generations, the number of migrants increases to Nm = 1, and the mean allele frequency is quickly pulled down towards that in the gene pool, to 0.155. The fixed-state approximation is the limit of low migration, and so is independent of Nm (appendix Ab figure 7). Indeed, allele frequencies are now often intermediate, and so this approximation fails (orange versus red, figure 2, middle). Nevertheless, it does give the important intuition that rates of change are proportional to migration, which is now m = 0.01, implying a 1/m ∼ 100 generation timescale for response of the population mean. In this model, variance is maintained by migration, and so the response to selection is proportional to m: we can see this in equation (2.3), where rates of change are proportional to m, for given Ns. After Nm returns to the original low value at 30 000 generations, there is a slow return to the original bimodal distribution, again captured by the fixed-state approximation (orange versus red at right of figure 2). Figure 2b shows the response to changes in Ns, which could arise through changes in selection strength, and/ royalsocietypublishing.org/journal/rstb Phil. Trans. R. Soc. B 377: 20210009 or changes in local effective population size. In this example, Nm = 0.05 throughout, and so the fixed-state approximation is accurate (orange versus red curves). The timescale is again set by m, which determines the rate at which variation is introduced into local demes. Since m = 5 × 10 −4 , it takes thousands of generations for the proportion of loci fixed for the '1' allele to respond to changes in selection strength. Figure 3 shows the time taken for a population to respond to changes in Nm (figure 3a) or Ns (figure 3b), as a function of the other parameter. In each figure, the two curves show the time to respond to changes in either direction. As we saw in figure 2a, an increase in Nm (lower plot of figure 3a) causes a much faster response than a decrease (upper plot of figure 3a), simply because high gene flow introduces more genetic variance. However, if selection is very strong, the response time becomes similar in either direction, and decreases in proportion to Ns (right of figure 3a). The response to changes in Ns take somewhat longer for an increase than a decrease (compare upper versus lower plot of figure 3b), but the main pattern here is that the response time decreases in proportion to Nm.

(b) Evolution of a metapopulation
We begin by considering the stationary state of a metapopulation, extending Szép et al. [8] by allowing a finite number of demes-in which case, a low rate of mutation is required to maintain variation in the long term. We then give an example that shows how variation is lost, as loci fix across the whole metapopulation. Finally  i þ Nmð1 À 2rÞ; the first term is derived from the fixed-state approximation, in the limit Nm ≪ 1, and the second from the deterministic model, which requires s > m(1 − 2ρ) for polymorphism (see [8, p. 1037 top left]). In a metapopulation with a finite number of demes, variation must ultimately be lost: we must include mutation to allow a non-trivial stationary state. In this section, we examine how the outcome depends on the relative rates of selection and drift (Ns) and on the relative rates of mutation and migration (μ/m). In particular, we show that with sufficiently many demes, the outcome is insensitive to the mutation rate. Figure 4 shows the stationary state in the limit of small Nm, derived using the fixed-state approximation. The top row of figure 4 shows how the fraction of demes fixed at equilibrium in a rare and common habitat (i.e. E½k 1 =n 1 and E½k 2 =n 2 , respectively) depend on the strength of selection when mutation is appreciable (figure 4a) compared to when it is weak ( figure 4b), where k i , i = 1, 2, are as defined earlier, n 1   and n 2 are the total number of demes in the rare and common habitats, respectively, n = n 1 + n 2 is the total number of demes in the metapopulation and mutation is assumed symmetric. Out of all demes, n, in the metapopulation, the focal allele is favoured in 20% of demes (constituting the rare habitat) so that n 1 = 0.2 n, and disfavoured twice as strongly in the remaining 80% of demes (constituting the common habitat) so that n 2 = 0.8 n. The blue (red) colour represents dynamics in the rare (common) habitat when we have a finite metapopulation (i.e. n = 50, 100, 200 and 400 demes represented by the solid, dashed, dotted and dot-dashed lines respectively). Whereas, the purple (orange) colour represent dynamics in the rare (common) habitat when we have an infinite metapopulation. When mutation is appreciable (μ/m = 0.05, figure 4a), the focal allele is unlikely to be lost by chance (blue colours). Furthermore, the equilibrium is insensitive to the number of demes, and close to the solution for an infinite population (as can be seen from the indistinguishability of the various lines for different N). When selection is strong (right of figure 4a), all demes are fixed for the favoured allele, whereas when selection is negligible, on average half of the demes are fixed for each allele. In between (i.e. 0.1 < Ns < 1), the allele favoured in the rare habitat becomes rare, being pulled to low frequency by migration from the commoner habitat, where it is more strongly disfavoured. When mutation is weak relative to migration (as is likely in nature, i.e. μ/m = 0.0005; figure 4b), this pattern is exaggerated. Above a critical value, Ns crit ∼ (1/ 2)log [(1 − ρ)/ρ] ∼ 0.7, polymorphism can be maintained by divergent selection, despite drift and gene flow. The equilibrium for an infinite population ( purple solid line, indicated by arrow) gives an upper bound, but stochastic loss from a finite set of demes reduces the expected frequency, and increases the critical Ns crit (solid, dashed, dotted and dotdashed blue lines around Ns ∼ 1, for 50, 100, 200 and 400 demes, respectively). There is a wide region (0.03 < Ns < 0.7) where the allele is almost absent, being swamped by gene flow. However, for very weak selection (figure 4b, l.h.s.), the frequency of the allele increases towards the symmetric neutral equilibrium, i.e. 0.5. In this regime, the frequencies in the two habitats are almost identical, and cannot be distinguished in the figure. Furthermore, in this regime (left of figure 4a), although selection is negligible within demes (Ns < 0.1), migration is much faster than mutation, and so selection over the whole metapopulation is effective in eliminating the allele that is deleterious in most demes. Therefore, in this scenario where mutation is rarer than migration, and selection is weak relative to drift within a single deme (i.e. μ ≪ m, Ns < 0.1), selection is nevertheless still effective at the level of the whole metapopulation and is especially so in the habitat which has more demes (left of figure 4a).
The bottom row of figure 4 shows the dependence on the relative rates of mutation versus migration, μ/m. With high mutation rates, the equilibrium approaches a fraction E½k=n ¼ 1=ð1 þ e À2Ns Þ, given by the ratio of fixation probabilities (equation (2.2)) in the fixed-state approximation. There is strong divergence when Ns 1 = 1 (right of figure 4c), and weaker divergence when selection is weak (figure 4d, Ns 1 = 0.1). With moderately strong selection (figure 4c), the allele that is less favoured overall is lost from the common habitat, independent of the number of demes and mutation rate (orange line indicated by arrow). In the rare habitat, with weak mutation (left of figure 4c), the locally favoured allele can be fixed in nearly half the demes in an infinite metapopulation (purple solid line indicated by arrow), but tends to be lost by chance from finite metapopulations, even with several hundred demes (solid, dashed and dotted blue lines). When selection is weak relative to local deme size (figure 4d ), selection can still be effective over the whole metapopulation, eliminating the allele that is disfavoured overall (left of figure 4d). However, when mutation becomes comparable with migration, polymorphism is maintained by mutation pressure, with some bias between habitats caused by weak selection (right of figure 4d).
We focus on the regime where selection is comparable to drift (Ns 1 ∼ 1), and mutation is weak. This corresponds to the right half of figure 4b (0.1 < Ns 1 ), and the middle of figure 4c (10 −4 < μ/m < 0.1). Then, as long as mutation is not extremely small, and there are enough demes, the stationary state is close to that in an infinite metapopulation (compare blue dashed and dotted lines with purple line in figure 4c). However, note that with weak mutation (μ/m ∼ 10 −4 − 10 −3 , say), the locally favoured allele tends to be lost even when there are several hundred demes.

(ii) Loss of diversity in a finite metapopulation
When deme sizes are fixed, and numbers of migrants are low enough that loci are typically fixed for one or other allele, the state of the metapopulation at each locus can be described by the number of demes, k i , in each habitat, i, that are fixed for the '1' allele. The distribution of k i evolves according to a transition matrix, and each locus follows an independent realization of the same stochastic process. Figures 8 and 9 in appendix Ac compare the dynamics of the fixed-state approximation with simulations, to illustrate its accuracy. For Nm = 0.05 (figure 8), there is reasonable agreement between simulations and the fixed-state approximation and for lower Nm (i.e. Nm = 0.01), agreement is even more close (figure 9). In both cases, variation is lost faster than predicted by the fixed-state approximation (compare red and black lines of figures 8a and 9a), because migration tends to swamp adaptive divergence. The timescale is inversely proportional to m, and is therefore slow. Here, we are focusing on the slow loss of adaptation through random drift in small populations; with higher migration rates, swamping by gene flow causes additional, faster, degradation.
Note that because the number of demes is limited, and because each deme flips between fixation for alternative alleles, there is substantial variability in average allele frequency between loci (grey lines). Therefore, adaptation is lost slightly faster in a finite than in an infinite metapopulation (compare black and magenta lines in figure 8a,b, which both derive from the fixed-state approximation). Nevertheless, the overall mean, averaged over 40 loci, changes smoothly and predictably (red curves in figures 8 and 9). We assume no mutation, and so all variation will inevitably be lost. However, because the total population is large (i.e. 100 demes × 50 individuals per deme = 5000 individuals), and because very low migration rate increases the effective size of the whole metapopulation, loss across the whole metapopulation is extremely slow: none of the 40 loci fix during the 10 4 generations shown in figures 8 and 9 (grey lines).
(iii) Response to changing conditions across the metapopulation Figure 5 shows some examples of the response to a change in conditions when parameters change gradually in all demes of the rare habitat. We focus on the rarer habitat since we are royalsocietypublishing.org/journal/rstb Phil. Trans. R. Soc. B 377: 20210009 mostly concerned with adaptation in the rare habitat and its degradation by gene flow. Figure 5a,d,g, respectively, show the consequence of a gradual increase in Ns 1 on the distribution of allele frequencies in the rare habitat, the equilibrium allele frequency, p (averaged across all demes in the rare habitat), and the half time, T 1/2 , taken to reach this equilibrium mean frequency. We begin initially with selection comparable to drift (Ns = 1) and with a fraction of demes fixed for the locally favoured allele, in the proportions predicted by an infinite metapopulation. After several thousand generations when the local allele has reached an equilibrium, we gradually increase Ns 1 until a new equilibrium is obtained and do this continuously until all loci in the rare habitat reach (near) perfect adaptation (figure 5a, grey lines). With selection comparable to drift, there is a substantial variation across loci with the 40 individual loci following markedly different paths and a chance loss at six of the loci (due to drift). As selection becomes stronger, the remaining polymorphic loci rise in frequency. Although there is still considerable variation in the rates of increase across loci (20 000-50 000 generations, figure 5a), the overall mean p is approached quite smoothly (red lines). These new equilibrium mean values have a positive dependence on Ns 1 (figure 5d) and it takes a shorter time to reach a new equilibrium with stronger selection (figure 5g). Figure 5b,e,h correspondingly shows a similar scenario as the above but with Ns 1 now changing in the opposite direction (i.e. from high to low value). Initially, with selection much stronger than drift, all 40 loci are at a considerably high frequency (near 1) with less variability amongst them (l.h.s of figure 5b, grey lines) so that the equilibrium mean allele frequency ∼1 (figure 5e, r.h.s.). Also, because of this strong selection, the new equilibrium is approached rather fast (figure 5h, r.h.s.). As Ns 1 declines, there is an apparent increase in the variability among loci so that with Ns 1 = 1 drift becomes sufficient enough to cause loss at three quarters of the loci (r.h.s. of figure 5b, grey lines). Because only 10/40 of the loci remain polymorphic, the overall mean is ∼0.1 (figure 5e, l.h.s.). Figure 5c,f,i shows the response (of the distribution of allele frequency, the equilibrium mean frequency p and the half time to p) to an increase in Nm with Ns 1 , Ns 2 = 10, −20 throughout. Since selection is generally strong relative to drift (Ns ≫ 1) and gene flow is initially low (Nm = 0.05, figure 5c, l.h.s.), the rare habitat is initially perfectly adapted with p 0:99 (figure 5f, l.h.s.). However, as Nm increases, gene flow gradually degrades adaptation as the focal allele is rapidly swamped until it is lost from the population after ∼13 000 generations (figure 5c, r.h.s.). Since no loci remain polymorphic at this point, p 0 (figure 5f, r.h.s.). The half time to reach the new equilibrium values depends non-monotonically on Nm.
In figure 5a-c, the dynamics are closely predicted by the small Nm limit (compare red and black curves) which is based simply on the rates of substitution in both directions.

Discussion
Our analysis uses simulation, the diffusion approximation and the 'fixed-state' approximation to understand how drift degrades adaptation in a finite metapopulation as well as how a finite metapopulation changes through time, as it responds to changes in both local and global conditions. The 'fixed-state' approximation applies either where variation is due to mutation (when it is plausible that Nμ < 1 20 000 40 000 60 000 time 0.5 royalsocietypublishing.org/journal/rstb Phil. Trans. R. Soc. B 377: 20210009 within local demes, or even at the level of the whole population), or when variation is maintained by divergent selection across the whole metapopulation, but migration is low relative to drift (Nm < 1).
Using the fixed-state approximation, it was shown that when selection is weaker than drift (i.e. Ns < 1), polymorphism can only be maintained for a very narrow range of habitat proportions (see appendix Aa figure 6). However with strong selection, this range becomes much wider. Also, using both simulations and the fixed state approximation, we showed that when conditions in a single deme of the metapopulation change, the population responds on a short timescale of order 1/m, simply because in the regime we study, local genetic variance is maintained by migration. Variation may be temporarily lost as local conditions change, but can quickly be recovered. On the other hand, when conditions change across the whole metapopulation, variation that was maintained by divergent selection can be permanently lost, and is only slowly recovered by mutation. Even under constant conditions, variation at a locus can be lost by chance, unless there are a very large number of demes.
To simplify our analysis, we assumed an island model, with a large number of spatially equivalent demes (i.e. soft selection). This is unlikely to be the case in nature, but may nevertheless capture the behaviour of spatially extended populations if there is long-range migration, which can introduce locally adaptive alleles from a distant habitat. It may be that a leptokurtic dispersal distribution can allow efficient adaptation, if locally favoured alleles are not swamped, and yet can be recovered by occasional long-range migration [14,15].
Our analysis can however be further extended to hard selection, by including explicit density regulation; Szép et al. [8] showed that one can still apply the diffusion approximation, provided that growth rates are not too high. With hard selection, substitution rates now depend on deme size through Ns, and through the number of immigrant alleles, mNp, mNq. This dependence can be approximated by assuming that the population size is determined by the genetic load. Sachdeva et al. [16] and Szép et al. [8] refer to this as the 'semi-deterministic' approximation which is accurate when demographic stochasticity is weak. One can apply the 'fixed-state' approximation by further assuming that there are enough loci that the mean load is proportional to the mean across loci of the number of demes fixed for one or the other allele. The transition matrix can then be calculated as before, but is now a function of the population sizes in the two habitats, {N 1 , N 2 } which both depend on the current state via the load. The key assumption here is that with enough loci, the population sizes change almost deterministically, following the distribution of states across loci. One complication with hard selection however is the existence of multiple stable equilibria: changing conditions would not just cause equilibria to shift but also changes the rates of transitions between equilibria.
A key assumption in our analysis is that selection is directional: in a given environment, alleles experience a fixed selection pressure, which tends to drive out variation. More often, selection may favour an intermediate optimum for a quantitative trait, i.e. stabilising selection, such that when the mean is well-adapted, alleles are close to neutral. Our modelling framework can describe this case, but it is much more complex, since many different allele combinations can achieve the same optimum. However, if selection on each allele is weak (Ns < 1), then the infinitesimal model [17] applies, and can also describe the population dynamics [18]. Local adaptation may be possible under higher migration rates in such a regime.
We have shown that it is very difficult for directional selection to maintain local adaptations when selection is weak, relative to deme size (Ns ∼ 1). Migration must also be weak if it is not to swamp adaptation, in which case alleles are typically near to fixation, limiting the genetic variance available for adaptation and preventing the recovery of variation lost to random drift. This contrasts with global adaptation: selection can be effective across the whole metapopulation, even if selection is weaker than drift within local demes (Ns < 1), provided that there is sufficient gene flow (Nm > 1) [11]. Thus, we expect local adaptation to depend on relatively more strongly selected alleles than global adaptations [19].
In this work, we have also introduced a novel approach to understanding the dynamical evolution of metapopulations. Although the full behaviour requires simulation, the diffusion approximation allows the stationary state to be calculated, and identifies the key dimensionless parameters. Moreover, when migration is rare, we can use a fixed-state approximation that connects population genetics with models of adaptive walks [9].
Our work suggests several open questions that invite theoretical study. First, although we show that local adaptation requires that directional selection be stronger than drift within demes, that may not be the case with stabilizing selection. Under the infinitesimal model [17], genetic variance may be due to weakly selected alleles (Ns < 1), and yet can still sustain adaptation of polygenic traits. Second, local adaptation may be greatly impeded by hard selection: if maladapted populations collapse, they cannot be the site of future adaptation. Testing the theory in nature is challenging, because it requires measurement of fitness as well as genetic data. Nevertheless, it might be possible to find how manipulation of local deme size and gene flow (or natural variation in these parameters) alters fitness. If we know of loci responsible for local adaptation (e.g. [20,21]), then the theory developed here can be applied more directly, though in practice it would need to be modified to account for actual spatial structure.
In conclusion, the methods introduced in this study allow us to explore how species adapt to diverse environments: we can find how organisms expand their range to a wider span of environments, through local adaptation that is sustained by variation maintained across the whole metapopulation. Many questions remain to be studied within this framework: for example, how populations adapt to large numbers of diverse habitats; whether population regulation (hard selection) leads to a feedback that impedes adaptation; and whether genetic variation can be better sustained under stabilizing selection towards a varying optimum, rather than the directional selection studied here. We believe that, despite the absence of explicit spatial structure, this approach will be a fruitful way to better understand what limits the range of environments that a species can occupy.
Data accessibility. All scripts used in this study are openly accessible through https://github.com/StochasticBiology/boolean-efflux.git. The data and mathematica codes are provided in the electronic supplementary material [22].
royalsocietypublishing.org/journal/rstb Phil. Trans. R. Soc. B 377: 20210009 Figure 6. Bounds on the proportions of habitat 1, ρ, between which polymorphism is possible, as a function of the strength of selection in that habitat, Ns 1 . For a given Ns 1 , ρ has to lie between the two curves for polymorphism to be maintained. The three sets of bounds correspond to Ns 2 / Ns 1 = 0.5, 1, 2 (black, blue and purple, respectively). These results apply in the limit of low migration, and soft selection. (Online version in colour.)